In [1]:
%matplotlib inline
%load_ext autoreload
%autoreload 2
In [2]:
from __future__ import division
from functions import *
from utils import *
from gevi_classes import *

/Users/GP1514/.pyenv/versions/anaconda3-2.5.0/lib/python3.5/site-packages/matplotlib/__init__.py:1350: UserWarning:  This call to matplotlib.use() has no effect
because the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)
In [3]:
%%javascript
IPython.OutputArea.auto_scroll_threshold = 1000;
In [4]:
# instantiate utility class
gr = Graph()

Load DATA

Data collected in 5 X 1min, with 1min pause in between 2921 images obtained in 60s

In [5]:
username = os.path.expanduser('~').split('/')[-1]
if username == "GP1514":
    print("At Imperial")
    mouseAPath = '/Volumes/DATA/GEVI/MouseA/'
    mouseBPath = '/Volumes/DATA/GEVI/MouseB/'
    mouseM1217FPath = '/Volumes/DATA/GEVI/M1217F/'
    mouseM1223MPath = '/Volumes/DATA/GEVI/M1223M/'
else:
    print("Using laptop")
    mouseAPath = '/Users/guillaume/Projects/GEVI-DATA/2014 Oct 27/'
    mouseBPath = '/Users/guillaume/Projects/GEVI-DATA/2014 Oct 22/'
    mouseCPath = '/Users/guillaume/Projects/GEVI-DATA/2014 Oct 28/'

At Imperial
In [6]:
discard = loadDic()
In [7]:
discard
Out[7]:
{'MouseA': {3: [0, 1, 3], 6: [1]},
 'MouseB': {3: [1, 4]},
 'MouseC': {},
 'MouseM1217F': {},
 'MouseM1223M': {}}
In [8]:
mouseA = Mouse('mouseA', mouseAPath, [3,4,5,6],discard['MouseA'] )
# mouseA = Mouse('mouseA', mouseAPath, [3],discard['MouseA'] )
mouseB = Mouse('mouseB', mouseBPath, [2,3,4,5],discard['MouseB'] )

mouseM1217F = Mouse('mouseM1217F', mouseM1217FPath, range(1,14),discard['MouseC'] )
mouseM1223M = Mouse('mouseM1223M', mouseM1223MPath, range(21,36),discard['MouseC'] )
In [9]:
mouseA.loadData()
mouseB.loadData()
mouseM1217F.loadData(split=True)
mouseM1223M.loadData(split=True)
Out[9]:
1
In [10]:
# mouseB.experiments[1].repeats[0].setInfo('End of transition, then discard')
# mouseB.experiments[2].repeats[0].setInfo('Transition to anesthesia')
# mouseA.experiments[0].repeats[2].setInfo('Keep')
# mouseA.experiments[2].repeats[1].setInfo('Transition to desynchronization: short periods of silence trigger recovery of hemodynamics')
# mouseA.experiments[3].repeats[1].setInfo('Woke up - Discard')
# mouseA.experiments[3].repeats[2].setInfo('Discard')
# mouseA.experiments[3].repeats[3].setInfo('Discard')
# mouseA.experiments[3].repeats[4].setInfo('Discard')
# mouseA.getInfo()
# mouseB.getInfo()
In [11]:
#len(mouseM1217F.experiments[0].repeats[0].mRatio)
len(mouseM1217F.experiments[0].repeats[0].ratio)
Out[11]:
11800

Plot Data

Mouse M1217F

In [12]:
gr.plotHV(mouseM1217F)

Mouse M1223M

In [13]:
gr.plotHV(mouseM1223M)

Transfer functions

$V * \alpha = H$

$\mathscr{F}$ transform : $\hat{V}.\hat{\alpha} = \hat{H}$

$\alpha = \mathscr{F}^{-1} \frac{\hat{H}}{\hat{V}}$

ALPHA CURVE FITTING

$\alpha(t)=\frac{\tau_0}{t+0.1} - \frac{t}{\tau_1}*e^{\tau_3-\frac{t}{\tau_2}}$

Extra functions

In [14]:
def meanCorH(mouse):
    corH = np.mean(flatten(mouse.experiments.corH))
    return corH

def meanCorR(mouse):
    corR = np.mean(flatten(mouse.experiments.corR))
    return corR
    
# optimization for mouseA
def fn(p_est):
    alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
    gr.computeCorr(mouseM1223M, alpha)
    corR = np.mean(flatten(mouseM1223M.experiments.corR))
    res = (1 - corR)**2
#     print('%.2f \t\t t0:%.5f \t t1:%.5f \t t2:%.5f \t t3:%.5f'%(res, p_est[0],p_est[1],p_est[2], p_est[3]))
    return res

def fnH(p_est):
    alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
    gr.computeCorr(mouseM1223M, alpha)
    corH = np.mean(flatten(mouseM1223M.experiments.corH))
    res = (1 - corH)**2
#     print('%.2f \t\t t0:%.5f \t t1:%.5f \t t2:%.5f \t t3:%.5f'%(res, p_est[0],p_est[1],p_est[2], p_est[3]))
    return res

def computeAndPlotCorr(mouse, alpha):
        for i,exp in enumerate(mouse.experiments):
            for j,rep in enumerate(exp.repeats):
                gr.plotRModel(mouse, exp = i, rep = j,alpha = alpha )
                gr.plotHModel(mouse, exp = i, rep = j,alpha = alpha )
                gr.corrCoeff(mouse, exp = i, rep = j,alpha = alpha )
                
x = np.real(xax(mouseM1217F.experiments[0].repeats[0].mRatio, 60000*4))
In [15]:
r= mouseA.experiments[0].repeats[0].ratio
out=np.ones((60,60))
In [ ]:
# filtering parameters
# mouseM1217F.minFreqAlpha = 1
# mouseM1217F.maxFreqAlpha = 200
# mouseM1223M.minFreqAlpha = 1
# mouseM1223M.maxFreqAlpha = 200
mouseA.window = 1
gr.fHmax = 30
gr.fRmax = 200

Tranfer functions mouse M1217F

In [33]:
gr.plotTF(mouseM1217F)
Out[33]:
1
In [35]:
gr.plotmTF(mouseM1217F)
Out[35]:
1

Transfer functions mouse M1223M

In [37]:
gr.plotTF(mouseM1223M)
Out[37]:
1
In [39]:
gr.plotmTF(mouseM1223M)
Out[39]:
1

Voltage and Hemo reconstruction : with mean alpha function

/!\ Models from the same mouse, different experiments: Model from M1217F, exp11, applied on M1217F, exp3

In [41]:
# alpha = mouseA.experiments[3].repeats[3].meanAlphaModel
# mouseM1217F.experiments[0].repeats[0].getAlpha()
alpha = mouseM1217F.experiments[10].repeats[0].alpha
gr.plotRModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha )
gr.plotHModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha )

(393.852692246+336.518394168j)

Model from M1223M, exp1, applied on M1217F, exp3

In [43]:
alpha = mouseM1223M.experiments[0].repeats[0].meanAlphaModel
gr.plotRModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha )
gr.plotHModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha )

(309.630249021+194.699778108j)

Voltage and Hemo reconstruction : with mean alpha model function

Model from M1223M, exp5, applied on M1217F, exp3

In [44]:
# alpha = mouseA.experiments[3].repeats[3].meanAlphaModel
# gr.fHmax=50
# gr.fRmax=50
# mouseM1223M.experiments[2].repeats[0].getAlpha()
alpha = mouseM1223M.experiments[4].repeats[0].meanAlphaModel
gr.plotRModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha, window=10 )
gr.plotHModel(mouseM1217F, exp = 2, rep = 0,alpha = alpha, window=10 )
plt.figure()
plt.plot(alpha)

(-197.579796829+63.6037191247j)
Out[44]:
[<matplotlib.lines.Line2D at 0x12532af98>]

Distribution of alpha model parameters $\tau_0, \tau_1, \tau_2, \tau_3$

$v * \alpha = h$

$\alpha(t)=\frac{\tau_0}{t+0.1} - \frac{t}{\tau_1}*e^{-\frac{t}{\tau_2}}$

In [45]:
alpha = mouseM1223M.experiments[2].meanAlphaModel
fontsize=16
plt.plot(xax(alpha,60000),alpha)
plt.xlabel('Time [s]')
plt.title('Alpha model of mean TF mouse B')
plt.text(30,-0.01, r'$ \tau_{0} = %.2g$' %mouseM1223M.experiments[2].meanAlphaParams[0], fontsize =fontsize)
plt.text(30,-0.02, r'$ \tau_{1} = %.2g$' %mouseM1223M.experiments[2].meanAlphaParams[1], fontsize =fontsize)
plt.text(30,-0.03, r'$ \tau_{2} = %.2g$' %mouseM1223M.experiments[2].meanAlphaParams[2], fontsize =fontsize)
Out[45]:
<matplotlib.text.Text at 0x1219f8400>
In [46]:
gr.computeCorr(mouseM1217F, alpha=None)
gr.computeCorr(mouseM1223M, alpha=None)
gr.computeAlphaModels(mouseM1217F)
gr.computeAlphaModels(mouseM1223M)
gr.plotParamsIndex(mouseM1217F)
gr.plotParamsIndex(mouseM1223M)
gr.plotParamsCor(mouseM1217F)
gr.plotParamsCor(mouseM1223M)
print([meanCorR(mouseM1217F),meanCorH(mouseM1217F),meanCorR(mouseM1223M), meanCorH(mouseM1223M)])

[0.1563448269735376, 0.26246088437116288, 0.027636849236815739, 0.31536924381042847]

Optimize for hemo

In [47]:
resultHemo = minimize(fnH, x0=[1,1,1], method = 'COBYLA',
                      options={'gtol': 1e-8, 'disp': True})
print(resultHemo)

     fun: 0.38823345491806766
   maxcv: 0.0
 message: 'Optimization terminated successfully.'
    nfev: 125
  status: 1
 success: True
       x: array([ 2.66382118,  2.27147196,  3.84925794])
In [48]:
p_est = resultHemo.x
alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
plt.plot(xax(alpha,4*60000),alpha)
# alpha2 = mouseB.experiments[2].meanAlphaModel
# plt.plot(xax(alpha2,60000),alpha2)
Out[48]:
[<matplotlib.lines.Line2D at 0x1209d8b38>]
In [49]:
# computeAndPlotCorr(mouseA, alpha)
gr.plotRModel(mouseM1217F,2,0,alpha)
gr.plotHModel(mouseM1217F,2,0,alpha)
# gr.plotHModel(mouseA,2,4,mouseB.experiments[2].meanAlphaModel)
# gr.plotRModel(mouseA,2,4,mouseB.experiments[2].meanAlphaModel)
def printl(tab):
    for t in tab:
        print(t)

(-292.841011211+187.05957221j)

Optimize for voltage

In [50]:
# bounds = [(-10,10),(-10,10),(-10,10),(-10,10)]
# resultVolt = differential_evolution(fn, bounds, init = 'random')
resultVolt = minimize(fn, method = 'COBYLA', x0=[-1e-2,1,10,1],
                      options={'gtol': 1e-6, 'disp': True})
print(resultVolt)

     fun: 0.64420198592358024
   maxcv: 0.0
 message: 'Optimization terminated successfully.'
    nfev: 602
  status: 1
 success: True
       x: array([  2.20903584,   1.60786837,  30.16489481,   0.6832862 ])
In [51]:
p_est = resultVolt.x
print(p_est)
alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]

[  2.20903584   1.60786837  30.16489481   0.6832862 ]
In [59]:
plt.title('Optimized Alpha model of mean TF mouse B, exp %d'%(2+mouseB.start))
plt.text(30,0.5, r'$ \tau_{0} = %.2g$' %p_est[0], fontsize =fontsize)
plt.text(30,1.5, r'$ \tau_{1} = %.2g$' %p_est[1], fontsize =fontsize)
plt.text(30,2.5, r'$ \tau_{2} = %.2g$' %p_est[2], fontsize =fontsize)
plt.plot(xax(alpha,4*60000),alpha)
Out[59]:
[<matplotlib.lines.Line2D at 0x156869198>]
In [60]:
gr.plotRModel(mouseM1217F,2,0,alpha)
gr.plotHModel(mouseM1217F,2,0,alpha)
##
gr.plotRModel(mouseM1217F,9,0,alpha)
gr.plotHModel(mouseM1217F,9,0,alpha)

(-181.051930801+54.0245725222j)
(408.654365153+18.9928711968j)

With the last model developed for mice A and B (from optimizing voltage model on voltage data from mouse B)

In [61]:
p_est = [0.4121683, 1.09384498, 3.68435018,  0.98029149]
In [62]:
alpha = [guess_function(xi,p_est[0],p_est[1],p_est[2]) for xi in x]
gr.plotRModel(mouseM1217F,2,0,alpha)
gr.plotHModel(mouseM1217F,2,0,alpha)

(757.701087231+84.1674460052j)
In [63]:
gr.plotRModel(mouseM1217F,9,0,alpha)
gr.plotHModel(mouseM1217F,9,0,alpha)

(460.342689682+70.7657955747j)

TF on different experiments: M1217, exp 10

Model from M1223M, exp1, applied on M1217F, exp10

In [64]:
alpha = mouseM1223M.experiments[0].repeats[0].meanAlphaModel
plt.figure(figsize=(10,1))
plt.plot(alpha)
gr.plotRModel(mouseM1217F, exp = 9, rep = 0,alpha = alpha )
gr.plotHModel(mouseM1217F, exp = 9, rep = 0,alpha = alpha )

(118.035268091-90.008853263j)
In [65]:
alpha = mouseM1223M.experiments[4].repeats[0].meanAlphaModel
plt.figure(figsize=(10,1))
plt.plot(alpha)
gr.plotRModel(mouseM1217F, exp = 9, rep = 0,alpha = alpha )
gr.plotHModel(mouseM1217F, exp = 9, rep = 0,alpha = alpha )

(408.654365153+18.9928711968j)
In [ ]: